Implicit solution algorithms are widely used to solve nonlinear time dependent transport problems in computational fluid dynamic (CFD) applications such as reservoir simulations. These implicit solution algorithms are commonly used to overcome time step size limitations of explicit methods. Often a solution to these transport problems for a new time level is obtained utilizing the well-known Newton-Raphson method.
FIG. 1 schematically illustrates a simple execution of the Newton-Raphson method for solving, or at least approximating, the solution for a function F=ƒ(S). FIG. 1 depicts a graph of the function:F=ƒ(S).  (1)
Solving ƒ(S)=F=0 determines a root, i.e., where the function ƒ(S) intersects with the horizontal axis, i.e., F=0. This intercept location for S is designated with the letter A. If this location, or a satisfactorily close approximation there to, can be determined, then Eqn. (1) is deemed to be solved.
The Newton-Raphson method requires making an initial guess for the solution. The initial guess is labeled as S0 in FIG. 1. A corresponding point B on the curve F=ƒ(S) has the coordinate (S0,ƒ(S0)). A tangent line is drawn at point B which intercepts with the horizontal axis. The slope of this tangent line is ƒ′(S0). Algebraically, the equation of the tangent at a general point S0 is:F−ƒ(S0)=ƒ′(S0)(S−S0)  (2)
The solution S is the horizontal coordinate of the tangent line when F=0. Solving for this coordinate:S=S0−ƒ(S0)/ƒ′(S0)  (3)
This coordinate is then used as a second guess or approximation for the solution of ƒ(S)=0 and is designated by S1. Accordingly,S1=S0−ƒ(S0)/ƒ′(S0)  (4)
A better estimate of the root of Eqn. (1) can be determined by again taking a tangent line at the coordinate (S1,ƒ(S1)), i.e., point C, and determining its intercept with the horizontal axis. This is accomplished using:S2=S1−ƒ(S1)/ƒ′(S1)  (5)
It can be seen that S2 is a better approximation than S1 for the solution of Eqn. (1). An even better estimate is made by taking the tangent to point D, (S2,ƒS2)), which results in the intercept at S3. In general, if an initial approximation or guess for the solution of F=ƒS) is Sn the next approximation Sn+1 may be determined using:Sn+1=Sn−ƒ(Sn)/ƒ′(Sn)  (6)
Eqn. (6) provides the basic formula for the Newton-Raphson method. By repeatedly applying the Newton-Raphson Eqn. (6) to a function ƒS), a better and better approximation for the solution to Eqn. (1) can be determined. However, this assumes that application of the Newton-Raphson method to find a solution on ƒS) is unconditionally stable. This is not always the case.
In reservoir simulation, the function to be solved is often a transport problem which is linearized around the best estimate of the solution and solved implicitly in each iteration. The transport problem is solved for each cell in a reservoir model in each Newton-Raphson iteration in a typical reservoir simulator. Once a converged solution is obtained after a number of Newton-Raphson iterations for a time step, the simulation then proceeds with the next time step. This approach is straightforward and unconditionally stable, if the flux function is convex (shock) as shown in FIG. 1, concave (rarefaction wave), or linear (contact discontinuity), which is the case for most computational fluid dynamics (CFD) problems. However, for S-shaped flux functions, which are encountered in most oil reservoir simulations, the Newton-Raphson method does not converge if the time steps selected are too large. In reservoir simulation, this problem is usually overcome by empirical time step control techniques. However, often the resulting time step size is conservative or too large, resulting in redundant computation or in time step cuts.
Accordingly, there is a need for an implicit solution algorithm for solving transport problems with S-shaped flux functions which is unconditionally stable, regardless of time step size. The present invention provides a solution to this need.